#clean up
rm(list=ls())

#packages
#library(XLConnect)
library(irr)
library(tools)
library(gmodels)
library(maps)
library(foreign)
library(lattice)
library(plm)
library(lmtest)
library(car)
library(xtable)
library(systemfit)

#load data
#setwd('/Volumes/MONOGAN/immigration/stateLaws/')
#setwd('/Users/jamie/Documents/immigration/stateLaws/')
data<-read.csv("immLawsMasterClean.csv")
state.aggregate<-read.csv("stateAggregate.csv")
policy.data<-read.csv("policyArea.csv")
panel.data<-read.csv("panelImmig.csv")

####INTERCODER RELIABILITY####
#coding years
table(data$year[!is.na(data$pro1)])
table(data$year[!is.na(data$pro2)])
table(data$year[!is.na(data$pro3)])
table(data$year[!is.na(data$pro4)])
table(data$year[!is.na(data$pro5)])
table(data$year[!is.na(data$pro6)])
table(data$year[!is.na(data$pro7)])

#How did all seven do in 2008?
data.scope.2008<-t(as.matrix(subset(data,select=c(scope1,scope2,scope3,scope4,scope5,scope6,scope7),subset=year==2008)))
kripp.alpha(data.scope.2008,"ordinal") #.715

data.pro.2008<-t(as.matrix(subset(data,select=c(pro1,pro2,pro3,pro4,pro5,pro6,pro7),subset=year==2008)))
kripp.alpha(data.pro.2008,"nominal")#.8

#How did everyone do in all times?
data.scope.subset<-t(as.matrix(subset(data,select=c(scope1,scope2,scope3,scope4,scope5,scope6,scope7))))
kripp.alpha(data.scope.subset,"ordinal")$value#.789

data.pro.subset<-t(as.matrix(subset(data,select=c(pro1,pro2,pro3,pro4,pro5,pro6,pro7))))
kripp.alpha(data.pro.subset,"nominal")$value#.789

####OVERALL SUMMARIES####
#overall number
#pdf("numberLaws.pdf")
plot(table(data$year),xlab="Year",ylab="Number of Laws",type="l")
points(x=2013,y=411)
text("411",x=2013,y=411,pos=4)
points(x=2005,y=39)
text("39",x=2005,y=39,pos=1)
#dev.off()

#topics
table(data$subject)
table(data$subject,data$year)
#pdf("topics.pdf")
topics<-100*(table(data$subject))/sum(table(data$subject))
names(topics)<-toTitleCase(names(table(data$subject)))
topics<-topics[order(topics)]
names(topics)[12]<-"Law Enforcement"
par(mar=c(5.1, 8 ,4.1 ,3.1))
barplot(topics,xlab='Percent of Laws',xlim=c(0,38),horiz=T,cex.names=.8,las=1)
text(x=topics,y=c(.75,1.75,3,4.25,5.5,6.65,7.8,9,10.25,11.5,12.65,13.75,15),labels=paste(round(topics,2)),pos=4)
box()
#dev.off()

#scope versus tone
CrossTable(round(data$scope),data$pro,prop.t=F,prop.r=F,prop.c=T)

####TEMPORAL GRAPH####
pro.year<-as.matrix(by(data$scope[data$pro==1],INDICES=data$year[data$pro==1],FUN=sum)); pro.year
con.year<-as.matrix(by(data$scope[data$pro==0],INDICES=data$year[data$pro==0],FUN=sum)); con.year
score.year<-log((pro.year+1)/(con.year+1)); score.year
#pdf("timeWelcoming.pdf")
plot(y=score.year,x=2005:2016,type='l',xlab="Year",ylab="Average Immigrant Welcomingness Score")
#dev.off()

###SCATTERPLOT OF ACTIVITY AND TONE###
state.aggregate$total<-as.vector(by(data$scope,INDICES=data$state,FUN=sum))
state.aggregate$abbr<-c('AL','AK','AZ','AR','CA','CO','CT','DE','FL','GA','HI','ID','IL','IN','IA','KS','KY','LA','ME','MD','MA','MI','MN','MS','MO','MT','NE','NV','NH','NJ','NM','NY','NC','ND','OH','OK','OR','PA','RI','SC','SD','TN','TX','UT','VT','VA','WA','WV','WI','WY')
#pdf("productionTone.pdf")
plot(y=state.aggregate$score,x=(state.aggregate$total),type='n',xlab="Weighted Number of Laws",ylab="Policy Tone",xlim=c(0,500))
text(y=state.aggregate$score,x=(state.aggregate$total),labels=state.aggregate$abbr)
abline(h=0,col='gray60')
abline(v=0,col='gray60')
#dev.off()

####MAP OF STATES####
ak<-state.aggregate$score[2]
hi<-state.aggregate$score[11]
score.48<-state.aggregate$score[-c(2,11)][c(1:20,20,21:48)] #Cut AK & HI, duplicate MI for both parts.
#7 levels
#Legend
#data(state, package = "datasets")

#List of State Map Objects
states.to.map<-c("alabama",
"arizona",
"arkansas",
"california",
"colorado",
"connecticut",
"delaware",
"florida",
"georgia",
"idaho",
"illinois",
"indiana",
"iowa",
"kansas",
"kentucky",
"louisiana",
"maine",
"maryland",
"massachusetts:main",
"michigan:north","michigan:south",
"minnesota",
"mississippi",
"missouri",
"montana",
"nebraska",
"nevada",
"new hampshire",
"new jersey",
"new mexico",
"new york:main",
"north carolina:main",
"north dakota",
"ohio",
"oklahoma",
"oregon",
"pennsylvania",
"rhode island",
"south carolina",
"south dakota",
"tennessee",
"texas",
"utah",
"vermont",
"virginia:main",
"washington:main",
"west virginia",
"wisconsin",
"wyoming")

#Policy Map
X<-score.48
splits<-quantile(X,c(0.1428571, 0.2857143, 0.4285714, 0.5714286, 0.7142857, 0.8571429))
n.col <- 7
br <- c(-999,splits, 999)
#shading<-c('dark red', 'red','pink','white','lightblue','blue', 'dark blue')
shading<-c('gray100', 'gray90', 'gray75', 'gray60', 'gray45', 'gray30', 'gray15')
X.grp<-findInterval(X, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
X.shad<-shading[X.grp]

#postscript("/home/monogan/DISSERTATION/document/leg_fig/policyMap2.eps", horizontal=FALSE, width=10, height=10, paper="special", onefile=FALSE, family='ComputerModern', pointsize=14)
#pdf("policyMap.pdf",width=10,height=10,pointsize=18)
#pdf("policyMapColor.pdf",width=10,height=10,pointsize=18)
map("state", fill=T, plot=T, col=X.shad, interior=T, region=states.to.map,ylim=c(25.12993,49.38323), xlim=c(-124.68134,-67.00742))#exact=T
box()
par(plt=c(.05,.35,.05,.25),new=TRUE)
plot.window(xlim=c(172.58641,229.98692),ylim=c(51.3522,71.29475))
map("world2",xlim=c(172.58641,229.98692),ylim=c(51.3522,71.29475),region=c("USA:alaska"), interior=T,fill=T,add=T,col=shading[findInterval(ak, vec=br, rightmost.closed = TRUE, all.inside = TRUE)],plot=T)#, exact=T)
box()
par(plt=c(.35,.60,.05,.25),new=TRUE)
plot.window(xlim=c(199.75,205.13580),ylim=c(19.00501,22.5))
map("world2",xlim=c(199.75,205.13580),ylim=c(19.00501,22.5),region=c("USA:Hawaii"), interior=T,fill=T,add=T,col=shading[findInterval(hi, vec=br, rightmost.closed = TRUE, all.inside = TRUE)],plot=T)#, exact=T)
box()
par(plt=c(.6,.95,.05,.25),new=TRUE)
plot.window(xlim=c(0,1),ylim=c(0,1))
legend("center", cex=.95, c('-1.38 - -0.46', '-0.46 - -0.06', '-0.06 -   0.19','  0.19 -   0.54','  0.54 -   0.91','  0.91 -   1.17','  1.17 -   2.13'),fill=shading,bty="n", ncol=1)
box()
#dev.off()

###PANEL MODEL###
panel.data$time<-panel.data$year-2005
panel.mod<-plm(score~cces.smooth+demUnif+repUnif+gsp+fb9015+fbPct2015+squireProfess+time,data=panel.data,model="random",index=c("state","year"));summary(panel.mod)
panel.t<-panel.mod$coefficients/sqrt(diag(panel.mod$vcov))
panel.mod.output<-cbind(panel.mod$coefficients,sqrt(diag(panel.mod$vcov)),panel.t,2*(1-pnorm(abs(panel.t))))
colnames(panel.mod.output)<-c("Coefficient","Std. Err.","t-ratio","p-value")
rownames(panel.mod.output)<-c("Intercept","Public Opinion Ideology","Unified Democrat","Unified Republican","Per Capita GSP","Delta Foreign Born 1990-2016","Foreign Born Percent 2015","Legislative Professionalism","Time")
xtable(panel.mod.output,digits=4,caption="Panel Model of Immigrant Policy Tone by State and Year, 2005-2015",label="panel.mod")

#some diagnostics
plot(y=resid(panel.mod),x=predict(panel.mod))
plot(y=resid(panel.mod),x=jitter(panel.data$time[!is.na(panel.data$score)]))
plot(y=resid(panel.mod),x=panel.data$cces.smooth[!is.na(panel.data$score)])
#alt<-plm(score~cces.smooth+demUnif+repUnif+gsp+fb9015+squireProfess+time,data=panel.data,model="random",subset=abs(cces.smooth)<50);summary(alt)
plot(y=resid(panel.mod),x=sqrt(panel.data$squireProfess[!is.na(panel.data$score)]))
plot(y=resid(panel.mod),x=panel.data$fb9015[!is.na(panel.data$score)])
bptest(panel.mod,studentize=F)
vif(panel.mod)

###SEPARATE PANEL MODELS###
panel.data.2<-subset(panel.data,subset=!is.na(pro) & !is.na(con))
panel.data.2$pro[panel.data.2$pro==0]<-0.1
panel.data.2$con[panel.data.2$con==0]<-0.1
r1<-log(pro)~cces.smooth+demUnif+repUnif+gsp+fb9015+fbPct2015+squireProfess+time
r2<-log(con)~cces.smooth+demUnif+repUnif+gsp+fb9015+fbPct2015+squireProfess+time
fitsur<-systemfit(list(pro.mod=r1,con.mod=r2),data=panel.data.2,method="SUR");summary(fitsur)
alt.pro<-plm(log(pro)~cces.smooth+demUnif+repUnif+gsp+fb9015+fbPct2015+squireProfess+time,data=panel.data.2,model="random",index=c("state","year"));summary(alt.pro)
alt.con<-plm(log(con)~cces.smooth+demUnif+repUnif+gsp+fb9015+fbPct2015+squireProfess+time,data=panel.data.2,model="random",index=c("state","year"));summary(alt.con)

z.1<-fitsur$coefficients[1:9]/sqrt(diag(fitsur$coefCov))[1:9]
z.2<-fitsur$coefficients[10:18]/sqrt(diag(fitsur$coefCov))[10:18]
p.1<-2*(1-pnorm(abs(z.1)))
p.2<-2*(1-pnorm(abs(z.2)))
sep.pan.out<-cbind(fitsur$coefficients[1:9],sqrt(diag(fitsur$coefCov))[1:9],p.1,fitsur$coefficients[10:18],sqrt(diag(fitsur$coefCov))[10:18],p.2)
colnames(sep.pan.out)<-c("Coef.","Std. Err.","p-value","Coef.","Std. Err.","p-value")
rownames(sep.pan.out)<-c("Intercept","Public Opinion Ideology","Unified Democrat","Unified Republican","Per Capita GSP","Delta Foreign Born 1990-2015","Foreign Born Percent 2015","Legislative Professionalism","Time")
xtable(sep.pan.out,digits=4,caption="Seemingly Unrelated Regression of Welcoming and Hostile Immigrant Policy Tone by State and Year, 2005-2016",label="sep.panel")


###CROSS-SECTIONAL MODEL###
cross.sec<-lm(score~cces.smooth+demUnif+repUnif+gsp+fb9015+fbPct2015+squireProfess,data=state.aggregate);summary(cross.sec)
xtable(cross.sec,digits=4,caption="Cross-Sectional Model of State Immigrant Policy Tone, Pooling All Laws from 2005-2015 Into One Measure",label="cross.sec")

#some diagnostics
plot(y=resid(cross.sec),x=predict(cross.sec))
plot(y=resid(cross.sec),x=state.aggregate$cces.smooth)
plot(y=resid(cross.sec),x=state.aggregate$squireProfess)
plot(y=resid(cross.sec),x=state.aggregate$fb9015)
plot(y=resid(cross.sec),x=state.aggregate$gsp)
bptest(cross.sec,studentize=F)
vif(cross.sec)


###SINGLE CODER MODEL###
table(data$pro6,useNA="always")
dim(state.aggregate)

pro.overall6<-as.matrix(by(data$scope6[data$pro6==1],INDICES=data$state[data$pro6==1],FUN=sum)); pro.overall6
con.overall6<-as.matrix(by(data$scope6[data$pro6==0],INDICES=data$state[data$pro6==0],FUN=sum)); con.overall6
state.names6<-row.names(pro.overall6)
score6<-log((pro.overall6+1)/(con.overall6+1))

state.aggregate6<-as.data.frame(cbind(state.names6))
state.aggregate6$score6<-log((pro.overall6+1)/(con.overall6+1))
names(state.aggregate6)[1]<-"state"
state.aggregate<-merge(x=state.aggregate,y=state.aggregate6,by="state")

single.mod<-lm(score6~cces.smooth+demUnif+repUnif+gsp+fb9015+fbPct2015+squireProfess,data=state.aggregate);summary(single.mod)
t1<-score~cces.smooth+demUnif+repUnif+gsp+fb9015+fbPct2015+squireProfess
t2<-score6~cces.smooth+demUnif+repUnif+gsp+fb9015+fbPct2015+squireProfess
fit.test<-systemfit(list(avg.mod=t1,single.mod=t2),data=state.aggregate,method="SUR");summary(fit.test)

r1<-c(1,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0)
r2<-c(0,1,0,0,0,0,0,0,0,-1,0,0,0,0,0,0)
r3<-c(0,0,1,0,0,0,0,0,0,0,-1,0,0,0,0,0)
r4<-c(0,0,0,1,0,0,0,0,0,0,0,-1,0,0,0,0)
r5<-c(0,0,0,0,1,0,0,0,0,0,0,0,-1,0,0,0)
r6<-c(0,0,0,0,0,1,0,0,0,0,0,0,0,-1,0,0)
r7<-c(0,0,0,0,0,0,1,0,0,0,0,0,0,0,-1,0)
r8<-c(0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,-1)
Rrestr<-rbind(r1,r2,r3,r4,r5,r6,r7,r8)
linearHypothesis(fit.test,Rrestr,test="F")


###SEPARATE AGGREGATE MODELS###
s1<-log(pro)~cces.smooth+demUnif+repUnif+gsp+fb9015+fbPct2015+squireProfess
s2<-log(con)~cces.smooth+demUnif+repUnif+gsp+fb9015+fbPct2015+squireProfess
fit.agg<-systemfit(list(pro.mod=s1,con.mod=s2),data=state.aggregate,method="SUR");summary(fit.agg)

z.1<-fit.agg$coefficients[1:8]/sqrt(diag(fit.agg$coefCov))[1:8]
z.2<-fit.agg$coefficients[9:16]/sqrt(diag(fit.agg$coefCov))[9:16]
p.1<-2*(1-pnorm(abs(z.1)))
p.2<-2*(1-pnorm(abs(z.2)))
sep.agg.out<-cbind(fit.agg$coefficients[1:8],sqrt(diag(fit.agg$coefCov))[1:8],p.1,fit.agg$coefficients[9:16],sqrt(diag(fit.agg$coefCov))[9:16],p.2)
colnames(sep.agg.out)<-c("Coef.","Std. Err.","p-value","Coef.","Std. Err.","p-value")
rownames(sep.agg.out)<-c("Intercept","Public Opinion Ideology","Unified Democrat","Unified Republican","Per Capita GSP","Delta Foreign Born 1990-2015","Foreign Born Percent 2015","Legislative Professionalism")
xtable(sep.agg.out,digits=4,caption="Seemingly Unrelated Regression of Welcoming and Hostile Immigrant Policy Tone, Pooling Across All Years for each State",label="sep.agg")

###GRAPH NUMBER OF WELCOMING AND HOSTILE LAWS BY PUBLIC IDEOLOGY###
pub.ideol<-seq(-34.59,15.03,by=.1)
pred.pro<-exp(cbind(1,pub.ideol,mean(state.aggregate$demUnif),mean(state.aggregate$repUnif),mean(state.aggregate$gsp),mean(state.aggregate$fb9015),mean(state.aggregate$fbPct2015),mean(state.aggregate$squireProfess))%*%fit.agg$coefficients[1:8])
pred.con<-exp(cbind(1,pub.ideol,mean(state.aggregate$demUnif),mean(state.aggregate$repUnif),mean(state.aggregate$gsp),mean(state.aggregate$fb9015),mean(state.aggregate$fbPct2015),mean(state.aggregate$squireProfess))%*%fit.agg$coefficients[9:16])

#pdf('ideolEffect.pdf')
plot(x=pub.ideol,y=pred.pro,type='l',ylim=c(0,70),lwd=2,xlab="Public Opinion Liberalism",ylab="Weighted Number of Laws")
lines(x=pub.ideol,y=pred.con,lty=2,col='red',lwd=2)
legend(x=-30,y=20,lty=c(1,2),col=c("black","red"),lwd=2,legend=c("Welcoming","Hostile"))
abline(h=0,col='gray60')
#dev.off()


###POLICY AREA MODELS###
table(policy.data$policy)
policy.wide<-reshape(policy.data,v.names=c("score","pro","con"),timevar="policy",idvar="state",direction="wide",sep=".")
a1<-score.appropriations~cces.smooth+demUnif+repUnif+gsp+fb9015+fbPct2015+squireProfess
a2<-score.benefits~cces.smooth+demUnif+repUnif+gsp+fb9015+fbPct2015+squireProfess
a3<-score.education~cces.smooth+demUnif+repUnif+gsp+fb9015+fbPct2015+squireProfess
a4<-score.employment~cces.smooth+demUnif+repUnif+gsp+fb9015+fbPct2015+squireProfess
a5<-score.health~cces.smooth+demUnif+repUnif+gsp+fb9015+fbPct2015+squireProfess
a6<-score.identification~cces.smooth+demUnif+repUnif+gsp+fb9015+fbPct2015+squireProfess
a7<-score.law~cces.smooth+demUnif+repUnif+gsp+fb9015+fbPct2015+squireProfess
a8<-score.licensing~cces.smooth+demUnif+repUnif+gsp+fb9015+fbPct2015+squireProfess
a9<-score.miscellaneous~cces.smooth+demUnif+repUnif+gsp+fb9015+fbPct2015+squireProfess
a10<-score.omnibus~cces.smooth+demUnif+repUnif+gsp+fb9015+fbPct2015+squireProfess
a11<-score.resolution~cces.smooth+demUnif+repUnif+gsp+fb9015+fbPct2015+squireProfess
a12<-score.trafficking~cces.smooth+demUnif+repUnif+gsp+fb9015+fbPct2015+squireProfess
a13<-score.voting~cces.smooth+demUnif+repUnif+gsp+fb9015+fbPct2015+squireProfess

policy.sur<-systemfit(list(mod.appropriations=a1,mod.benefits=a2,mod.education=a3,mod.employment=a4,mod.health=a5,mod.identification=a6,mod.law=a7,mod.licensing=a8,mod.miscellaneous=a9,mod.omnibus=a10,mod.resolution=a11,mod.trafficking=a12,mod.voting=a13),data=policy.wide,method="SUR");summary(policy.sur)

sink("policySUR.tex")
xtable(summary(policy.sur$eq[[1]])$coefficients,digits=4,caption="Model of Appropriations Policy in Seemingly Unrelated Regression",align="lrrrr")
xtable(summary(policy.sur$eq[[2]])$coefficients,digits=4,caption="Model of Benefits Policy in Seemingly Unrelated Regression",align="lrrrr")
xtable(summary(policy.sur$eq[[3]])$coefficients,digits=4,caption="Model of Education Policy in Seemingly Unrelated Regression",align="lrrrr")
xtable(summary(policy.sur$eq[[4]])$coefficients,digits=4,caption="Model of Employment Policy in Seemingly Unrelated Regression",align="lrrrr")
xtable(summary(policy.sur$eq[[5]])$coefficients,digits=4,caption="Model of Health Policy in Seemingly Unrelated Regression",align="lrrrr")
xtable(summary(policy.sur$eq[[6]])$coefficients,digits=4,caption="Model of Identification Policy in Seemingly Unrelated Regression",align="lrrrr")
xtable(summary(policy.sur$eq[[7]])$coefficients,digits=4,caption="Model of Law Enforcement Policy in Seemingly Unrelated Regression",align="lrrrr")
xtable(summary(policy.sur$eq[[8]])$coefficients,digits=4,caption="Model of Licensing Policy in Seemingly Unrelated Regression",align="lrrrr")
xtable(summary(policy.sur$eq[[9]])$coefficients,digits=4,caption="Model of Miscellaneous Policy in Seemingly Unrelated Regression",align="lrrrr")
xtable(summary(policy.sur$eq[[10]])$coefficients,digits=4,caption="Model of Omnibus Policy in Seemingly Unrelated Regression",align="lrrrr")
xtable(summary(policy.sur$eq[[11]])$coefficients,digits=4,caption="Model of Resolutions in Seemingly Unrelated Regression",align="lrrrr")
xtable(summary(policy.sur$eq[[12]])$coefficients,digits=4,caption="Model of Trafficking Policy in Seemingly Unrelated Regression",align="lrrrr")
xtable(summary(policy.sur$eq[[13]])$coefficients,digits=4,caption="Model of Voting Policy in Seemingly Unrelated Regression",align="lrrrr")
sink()
